สถิติแบบ Frequentist

สถิติแบบดั้งเดิมหรือที่เรียกกันว่า frequentist เป็นสถิติวิเคราะห์ที่ถูกใช้งานเป็นหลักในปัจจุบัน สถิติเกือบทุกตัวที่ผู้เรียนได้เรียนมาจนถึงปัจจุบันนี้ล้วนเป็นสถิติวิเคราะห์ภายใต้กระบวนทัศน์แบบ frequentist ทั้งสิ้น สถิติแบบ frequentist นี้มีข้อตกลงเบื้องต้นว่า (1) พารามิเตอร์เป็นค่าคงที่ ที่มีอยู่จริงในประชากรเป้าหมาย อย่างไรก็ตามด้วยข้อจำกัดที่ไม่สามารถเก็บรวบรวมข้อมูลจากทั้งประชากรได้ทำให้ผู้วิเคราะห์ไม่ทราบค่าพารามิเตอร์ดังกล่าว และจะประมาณโดยใช้ข้อมูลจากตัวอย่าง (2) ข้อมูลตัวอย่าง (sample data) เขียนแทนด้วย \(\bf{x}= \left \{\bf{x_1}, \bf{x_2},...,\bf{x_n} \right \}\) เป็นค่าที่สุ่มมาจากการประชากร ดังนั้นตัวอย่างที่สุ่มได้จึงเป็นค่าสุ่มที่และมีความเป็นไปได้ทั้งหมดเท่ากับ เมื่อ \(N\) คือขนาดประชากร และ \(n\) คือขนาดตัวอย่าง จากข้อตกลงเบื้องต้นข้อที่ (2) นี้จึงทำให้ (3) ค่าสถิติที่ประมาณได้จากข้อมูลตัวอย่างเป็นตัวแปรสุ่มที่มีการแจกแจงความน่าจะเป็น เรียกการแจกแจงความน่าจะเป็นของค่าสถิตินี้ว่า การแจกแจงความน่าจะเป็นของฟังก์ชันที่ได้จากตัวอย่างสุ่ม (sampling distribution)

ลักษณะการแจกแจงความน่าจะเป็นของค่าสถิติดังกล่าวมักประมาณได้ด้วย 2 วิธีการ วิธีการแรก ประมาณโดยใช้ทฤษฎีลิมิตลู่เข้าสู่ส่วนกลาง (central limit theorem) ซึ่งเป็นการคาดการณ์ลักษณะการแจกแจงของค่าสถิติในเชิงทฤษฎี เช่น การประมาณช่วงความเชื่อมั่นของค่าเฉลี่ยมีสูตรเป็น \(\overline{X}\pm t_{\alpha/2,n-1} \frac{S}{\sqrt{n}}\) จะเห็นว่าการประมาณช่วงความเชื่อมั่นดังกล่าวมีการใช้การแจกแจงที (student t distribution) ในการคำนวณส่วน margin of error การแจกแจงดังกล่าวพิสูจน์ในทางทฤษฎีโดยใช้ central limit theorem ได้ว่า เป็น sampling distribution ของ \(\overline{X}\) หากข้อมูลตัวอย่างที่ใช้ในการวิเคราะห์สุ่มมาจากประชากรที่มีการแจกแจงแบบปกติหรือใกล้เคียง แต่ไม่ทราบค่าความแปรปรวนของประชากร หรือถ้าหากประชากรไม่ได้มีการแจกแจงแบบปกติ แต่ตัวอย่างมีขนาดใหญ่เพียงพอ ก็ยังสามารถใช้ sampling distribution ดังกล่าวประมาณการแจกแจงของ \(\overline{X}\) ได้อยู่ หรือในการทดสอบสมมุติฐาน \(H_0: \mu=\mu_0\) สถิติทดสอบที่ใช้คือ \(t^*=\frac{\overline{X}-\mu_0}{S/ \sqrt{n}}\) จะมี sampling distributionในทางทฤษฎีเป็นการแจกแจงแบบที ที่มีองศาความเป็นอิสระเท่ากับ \(n-1\) เช่นเดียวกันหากข้อกำหนดเบื้องต้นของ central limit theorem เป็นจริงหรือใกล้เคียง การทราบการแจกแจงของค่าสถิติดังกล่าวทำให้ผู้วิเคราะห์สามารถประเมินความคลาดเคลื่อนในการตัดสินใจของการทดสอบสมมุติฐาน และสรุปผลได้ดังที่ได้เรียนไปแล้วในรายวิชาพื้นฐาน

ความถูกต้องของทฤษฎีนี้ขึ้นอยู่กับปัจจัยแวดล้อมที่เกี่ยวข้องว่ามีความสอดคล้องกับทฤษฎีมากน้อยเพียงใด ปัจจัยหนึ่งที่มีความสำคัญคือขนาดตัวอย่าง กล่าวคือหากขนาดตัวอย่างมีค่าน้อยเกินไป หรือการแจกแจงของข้อมูลมีความเบี่ยงเบนออกไปจากข้อตกลงเบื้องต้นของการวิเคราะห์ไปมาก ค่าสถิติที่คำนวณได้อาจไม่มีลักษณะการแจกแจงที่เป็นไปตามทฤษฎี และส่งผลให้การอนุมานเชิงสถิติมีความคลาดเคลื่อน อีกวิธีการหนึ่งที่สามารถใช้ในการประมาณ sampling distribution ของค่าสถิติได้คือใช้วิธีการสุ่มซ้ำ (resampling method) เช่น boostraping หรือ jackknife เพื่อประมาณแนวโน้มการแจกแจงของค่าสถิติ เป็นต้น

t.test(dat$Ach, mu=5.5)

## 
##  One Sample t-test
## 
## data:  dat$Ach
## t = 5.078, df = 149, p-value = 1.123e-06
## alternative hypothesis: true mean is not equal to 5.5
## 95 percent confidence interval:
##  5.709732 5.976934
## sample estimates:
## mean of x 
##  5.843333

คำถาม

1. output ข้างต้นแสดงผลการประมาณช่วงความเชื่อมั่น และทดสอบสมมุติฐานทางสถิติเกี่ยวกับผลสัมฤทธิ์ทางการเรียนของนักเรียน (คะแนนเต็ม 10) โดยจากการประมาณด้วยช่วงความเชื่อมั่น 95% ในข้างต้นพบว่ามีค่าเท่ากับ \([5.71, 5.98]\) ซึ่งหมายความว่าอะไร?

2. จากผลการทดสอบสมมุติฐาน \(H_0: \mu_{Ach}=5.5\) พบว่ามีค่า p-value < .000 หมายความว่าอะไร?

3. เราสามารถตัดสินใจยอมรับสมมุติฐาน \(H_0\) ได้หรือไม่?

จากข้อตกลงเบื้องต้นของสถิติแบบ frequentist นี้ทำให้การวิเคราะห์ข้อมูลเกิดข้อจำกัดต่าง ๆ ดังนี้

1. ด้านที่สำคัญคือด้านการแปลผลการวิเคราะห์ที่จะเห็นว่าความหมายของผลการวิเคราะห์ต่าง ๆ เป็นการกล่าวถึงพารามิเตอร์ที่สนใจแบบอ้อม ๆ ทั้งสิ้น ไม่มีการวิเคราะห์ใดที่สามารถอ้างอิงไปยังพารามิเตอร์ที่สนใจได้โดยตรง ทั้งนี้เป็นเพราะการอนุมานเชิงสถิติแบบ frequentist พึ่งพาเครื่องมือหลักคือ samping distribution ของค่าสถิติ ในขณะที่พารามิเตอร์ถูกสมมุติให้เป็นค่าคงที่ตั้งแต่ต้น probability statement ต่าง ๆ จึงเป็นของค่าสถิติทั้งสิ้น ดังความความเชื่อมั่นที่รับประกันในช่วงความเชือมั่น และค่าความคลาดเคลื่อนในการทดสอบสมมุติฐานจึงเป็นค่าความน่าจะเป็นที่ใช้รับประกันความเป็นไปได้ของค่าสถิติทั้งสิ้นไม่ใช่ของค่าพารามิเตอร์ที่สนใจจริง ๆ ในการวิเคราะห์

2. ด้านที่สองคือด้านการประมาณค่าพารามิเตอร์ภายใต้โมเดลซับซ้อน เครื่องมือสำคัญสำหรับประมาณค่าพารามิเตอร์ในโมเดลต่าง ๆ ของสถิติแบบ frequentist คือตัวประมาณ เช่น ตัวประมาณกำลังสองน้อยสุด ตัวประมาณภาวะความควรจะเป็นสูงสุด ซึ่งล้วนเป็นวิธีการที่จะหาค่าประมาณที่ดีที่สุดโดยอิงจากจุดต่ำสุด หรือจุดสูงสุดของฟังก์ชันวัตถุประสงค์ ดังตัวอย่างในรูปด้านล่าง ซึ่งในกรณีที่โมเดลมีความซับซ้อนมาก ๆ ฟังก์ชันวัตถุประสงค์ดังกล่าวจะมีความซับซ้อนโดยมีจำนวนมิติที่สูงขึ้น และอาจมีจุดสูงต่ำ และ/หรือ จุดต่ำสุดสัมพัทธ์จำนวนมาก ซึ่งเป็นอุปสรรคให้การประมาณค่าด้วยวิธีการในข้างต้นอาจได้ค่าประมาณที่คลาดเคลื่อนไปจากค่าที่เหมาะสมที่สุด หรืออัลกอริทึมไม่สามารถหาค่าประมาณที่เหมาะสมได้ภายใต้จำนวนรอบของการประมาณที่กำหนด

3. สืบเนื่องจากที่การอนุมานเชิงสถิติแบบ frequentist นั้นมักอิงกับทฤษฎีความน่าจะเป็นที่อิงข้อกำหนดเบื้องต้นที่ขนาดตัวอย่างต้องมีขนาดใหญ่เพียงพอ ซึ่งทำให้สถิติแบบ frequentist มักทำงานได้ไม่ดีในกรณีที่ตัวอย่างมีขนาดเล็ก

สถิติแบบเบส์ (Bayesian Statistics)

สถิติแบบเบส์เป็นกระบวนทัศน์ในการวิเคราะห์ข้อมูลอีกลักษณะหนึ่งที่มีความแตกต่างจากไปกระบวนทัศน์แบบดั้งเดิมที่นักวิเคราะห์ข้อมูลทั่วไปคุ้นเคยกัน ซึ่งทำให้สถิติแบบเบส์มีจุดเด่นและสามารถแก้ไขข้อจำกัดที่พบในสถิติแบบ frequentist ได้ ความแตกต่างดังกล่าวเริ่มตั้งแต่ข้อตกลงเบื้องต้นของสถิติแบบเบส์ที่แตกต่างจากสถิติแบบดั้งเดิมอย่างมาก กล่าวคือ

(1) พารามิเตอร์เป็นตัวแปรสุ่ม สถิติแบบเบส์มองว่าพารามิเตอร์เป็นค่าที่ผู้วิเคราะห์สนใจจะศึกษาและต้องการทราบอย่างแท้จริง แต่อย่างไรก็ตามเนื่องจากผู้วิเคราะห์ไม่ทราบข้อมูลประชากร การคาดการณ์หรือความเชื่อ (belief) เกี่ยวกับพารามิเตอร์ดังกล่าวจึงมีความไม่แน่นอน สภาพดังกล่าวจึงสมเหตุสมผลที่จะกำหนด probability statement ให้กับพารามิเตอร์ในโมเดล

จากข้อกำหนดนี้ผลการประมาณค่าพารามิเตอร์ที่ได้จากวิธีการแบบเบส์จึงจะไม่ได้ค่าประมาณแบบค่าเดียว แต่จะได้เป็นการประมาณการแจกแจงของพารามิเตอร์แทน เช่น การประมาณค่าเฉลี่ยผลสัมฤทธิ์ของนักเรียน ผลการวิเคราะห์แบบเบส์จะไม่ได้ให้เป็นค่าประมาณเพียงค่าเดียวคือ 5.84 คะแนน ดัง output ข้างต้น แต่จะให้เป็นการแจกแจงของพารามิเตอร์ค่าเฉลี่ยผลสัมฤทธิ์ทางการเรียน ซึ่งเรียกว่าการแจกแจงความน่าจะเป็นภายหลัง (posterior distribution) และสารสนเทศเกี่ยวกับค่าเฉลี่ยของผลสัมฤทธิ์ดังกล่าวจะถูกสกัดออกมาจากการแจกแจงนี้

(2) ข้อมูลตัวอย่างเป็นค่าคงที่ ถึงแม้ว่าในความเป็นจริงตัวอย่างที่สุ่มมาจากประชากรจะมีความเป็นไปได้ที่หลากหลาย และมีความไม่แน่นอนในการได้มาซึ่งตัวอย่างแต่ละชุด อย่างไรก็ตามเมื่อดำเนินการสุ่มตัวอย่างมาทำการวิเคราะห์แล้ว ความเป็นไปได้จำนวนมากก่อนการสุ่มตัวอย่างจะเหลือเพียงความเป็นไปได้เดียว สถิติแบบเบส์จึงถือว่าข้อมูลตัวอย่างเป็นค่าคงที่ ไม่มีการกำหนด probability statetment ให้กับข้อมูลตัวอย่างนี้

สถิติแบบเบส์สามารถนำมาใช้เพื่อแก้ไขข้อจำกัดที่เกิดขึ้นจากสถิติแบบ frequentist ได้ดังนี้

1. จากข้อกำหนดเบื้องต้นข้างต้น การอนุมานเชิงสถิติด้วยสถิติแบบเบส์จะดำเนินการผ่านการแจกแจงความน่าจะเป็นของพารามิเตอร์โดยตรงที่เรียกว่า การแจกแจงความน่าจะเป็นภายหลัง (posterior distribution) การแจกแจงดังกล่าวเป็นเครื่องมือสำคัญในการอนุมานเชิงสถิติแบบเบส์ แทน sampling distribution ในการอนุมานเชิงสถิติแบบดั้งเดิม ซึ่งจุดนี้ทำให้สถิติแบบเบส์สร้างความแตกต่างอย่างมากเมื่อเปรียบเทียบกันสถิติแบบดั้งเดิม โดยจะเห็นว่าการประมาณค่าพารามิเตอร์แบบเบส์ไม่ได้อิงกับค่าสถิติเพียงค่าเดียวเหมือนกับสถิติแบบดั้งเดิม แต่อิงกับการแจกแจงความน่าจะเป็นภายหลัง

2. การที่สถิติแบบเบส์ใช้การแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์เป็นเครื่องมือหลักในการอนุมานเชิงสถิติ จึงทำให้การวิเคราะห์ทำได้อย่างตรงไปตรงมาโดยไม่ต้องอาศัยทฤษฎีหรือวิธีการที่ซับซ้อนเหมือนกับสถิติแบบดั้งเดิม กล่าวคือผู้วิเคราะห์ใช้เพียงสถิติพื้นฐาน เช่น ค่าเฉลี่ย มัธยฐาน ฐานนิยม ส่วนเบี่ยงเบนมาตรฐาน หรือเปอร์เซ็นไทล์ เป็นเครื่องมือในการสรุปสารสนเทศเกี่ยวกับพารามิเตอร์ที่สนใจออกมาจากการแจกแจงความน่าจะเป็นภายหลังที่ประมาณได้ ก็เพียงพอแล้ว นอกจากนี้การแปลผลเพื่ออนุมานเกี่ยวกับพารามิเตอร์ยังทำได้ง่ายและตรงไปตรงมา และได้ข้อสรุปโดยตรงไปยังปริภูมิของพารามิเตอร์ที่เป็นเป้าหมาย แตกต่างจากข้อสรุปที่ได้จากสถิติแบบดั้งเดิมที่เป็นข้อสรุปเกี่ยวกับการแจกแจงของค่าสถิติซึ่งอยู่บนปริภูมิตัวอย่างเท่านั้น

3. การประมาณการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ในโมเดล จะใช้สารสนเทศสองส่วนด้วยกัน ส่วนแรกเป็นสมมุติฐานหรือความเชื่อเบื้องต้น (prior belief) เกี่ยวกับพารามิเตอร์ที่สนใจ ซึ่งกำหนดอยู่ในรูปการแจกแจงความน่าจะเป็นที่เรียกว่า การแจกแจงความน่าจะเป็นก่อนหน้าของพารามิเตอร์ (prior distribution) และส่วนที่สองคือสารสนเทศที่ได้จากข้อมูลตัวอย่าง ซึ่งแตกต่างจากสถิติแบบดั้งเดิมที่จะใช้สารสนเทศจากข้อมูลตัวอย่างเท่านั้น

ลองพิจารณาตัวอย่างต่อไปนี้

สมมุติว่าต้องการคาดการณ์คะแนนสอบของนายบูม ที่มีความเป็นไปได้ 4 แบบดังรูป

ก่อนการเก็บรวบรวมข้อมูลจริงในชั้นเรียนของนายบูม สมมุติว่าผู้วิเคราะห์มีข้อมูลในอดีตที่เกี่ยวข้องกับคะแนนสอบของนายบูมดังนี้

1. ไม่มีข้อมูลใด ๆ ที่เป็นประโยชน์ (ให้น้ำหนักกับความเป็นไปได้ A1, A2, A3 และ A4 อย่างไร?)

2. หากทราบเพิ่มเติมว่าการสอบที่ผ่าน ๆ มาโดยเฉลี่ยนายบูมสอบได้คะแนนคิดเป็นอันดับที่ 4 ของชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 จะเหมือนเดิมมั้ย?)

3. หากทราบเพิ่มเติมอีกว่า การสอบที่ผ่าน ๆ มา แก้วซึ่งเป็นเพื่อนในห้องเรียนของนายบูม สอบได้คะแนน 38 คะแนน คิดเป็นคะแนนที่น้อยที่สุดในชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 เหมือนเดิมหรือไม่?)

4. หากทราบอีกว่า ในการสอบที่ผ่าน ๆ มา นิดสอบได้คะแนนโดยเฉลี่ยคิดเป็น 89 คะแนน ซึ่งเป็นอันดับที่ 3 ของชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 เหมือนเดิมหรือไม่?)

การแจกแจงก่อนหน้าของคะแนนสอบนายบูม

ผู้วิเคราะห์สามารถปรับสารสนเทศในการแจกแจงความน่าจะเป็นก่อนหน้าในข้างต้นได้ หากมีการเก็บรวบรวมข้อมูลตัวอย่างเพิ่มเติม เช่น

สมมุติว่าการแจกแจงความน่าจะเป็นก่อนหน้าของคะแนนนายบูมเป็นแบบ uniform และมีข้อมูลการสอบครั้งล่าสุดเพิ่มเติมดังรูป

อีกตัวอย่างหนึ่ง สมมุติว่า การแจกแจงความน่าจะเป็นก่อนหน้าของคะแนนนายบูมมีน้ำหนักสูงที่ A3 และมีข้อมูลการสอบครั้งล่าสุดเพิ่มเติมดังรูป

กำหนดตัวแปรในการวิเคราะห์

สมมุติว่าผู้วิจัยต้องการวิเคราะห์ความลำเอียงในการออกหน้าหัวหรือก้อยของเหรียญที่ผลิตขึ้นมารุ่นหนึ่ง การวิเคราะห์นี้ตัวแปรที่สนใจจึงเป็นผลลัพธ์ที่สังเกตได้จากการโยนเหรียญในแต่ละครั้ง กำหนดให้ \(y_i\) แทนค่าสังเกตที่ได้จากการโยนเหรียญในครั้งที่ \(i\) จะได้ว่า

\(y_i=\left\{\begin{matrix} 1 & ; H\\ 0 & ; T \end{matrix}\right.\)

กำหนดโมเดลของค่าสังเกต

เมื่อกำหนดตัวแปรที่สนใจแล้วจะพบว่าตัวแปรที่สนใจมีเพียงตัวแปรตาม 1 ตัว โดยมีค่าสังเกตที่เป็นไปได้เพียง 2 ค่าได้แก่ H หรือ T ซึ่งแทนด้วย 1 กับ 0 ตามลำดับ โมเดลของค่าสังเกต \(y_i\) นี้จึงสามารถกำหนดได้ด้วยการแจกแจงความน่าจะเป็นแบบ Bernoulli กล่าวคือ ความน่าจะเป็นที่จะออกหน้าหัวมีค่าเท่ากับ \(p(y_i=1|\theta)=\theta\) และความน่าจะเป็นที่ออกหน้าก้อยมีค่าเท่ากับ \(p(y_i=0|\theta)=1-\theta\) ซึ่งเขียนเป็นสมการรวมได้ดังนี้

\(p(y_i|\theta)=\theta^{y_i}(1-\theta)^{1-y_i}\) เมื่อ \(i = 1,2,3,...,n\)

จากโมเดลในข้างต้นจะเห็นว่าการเกิดค่าสังเกต \(y_i\) ขึ้นอยู่กับพารามิเตอร์ \(\theta\) ที่มีค่าอยู่บนช่วง \([0,1]\) หากพารามิเตอร์ดังกล่าวมีค่าเท่ากับ \(\theta=0.5\) นั่นหมายถึงโอกาสในการออกหน้าหัวและก้อยมีค่าเท่ากัน แต่ถ้าหากพารามิเตอร์ \(\theta>0.5\) นั่นหมายถึงโอกาสในการออกหน้าหัวมีมากกว่า จากความหมายดังกล่าวทำให้สามารถมองได้ว่าพารามิเตอร์ \(\theta\) เป็นพารามิเตอร์ความลำเอียงของเหรียญ หากประมาณพารามิเตอร์ดังกล่าวได้จะสามารถสรุปความลำเอียงของเหรียญได้

เนื่องจากเราไม่ได้เก็บรวบรวมข้อมูลเพียงค่าสังเกตเดียว แต่เราทำการทดลองซ้ำ ๆ จำนวน \(n\) ครั้ง เนื่องจากการทดลองแต่ละครั้งเป็นอิสระซึ่งกันและกัน ทำให้โมเดลของชุดค่าสังเกตหรือเรียกอย่างเป็นทางการว่า ฟังก์ชันภาวะความควรจะเป็นดังนี้

\(p(\bf{y}|\theta)=p(y_1|\theta)p(y_2|\theta)...p(y_n|\theta)\)

ซึ่งมีค่าเท่ากับ

\(p(\bf{y}|\theta)=\theta^{\sum_{i=1}^ny_i}(1-\theta)^{n-\sum_{i=1}^ny_i}\)

ในฟังก์ชันภาวะความควรจะเป็นนี้จะให้ค่าความน่าจะเป็นหรือค่าความหนาแน่นของข้อมูลตัวอย่าง \(\bf{y}\) บนแต่ละค่าที่เป็นไปได้ของพารามิเตอร์ \(\theta\) เนื่องจากข้อมูลตัวอย่างเป็นค่าคงที่ ดังนั้นค่าความหนาแน่นนี้จึงเปลี่ยนแปลงไปตามค่าของพารามิเตอร์ \(\theta\) โดยหากค่าความหนาแน่นนี้มีค่าสูง ณ ค่า \(\theta\) ค่าใด นั่นหมายถึงว่าพารามิเตอร์ค่านั้นทำให้โมเดลของค่าสังเกตกับค่าสังเกตจริงมีความสอดคล้องกันสูง ข้อสังเกตนี้ถูกนำไปใช้เพื่อหาค่าประมาณด้วยวิธีภาวะความควรจะเป็นสูงสุด (maximum likelihood: ML) แต่ในวิธีการแบบเบส์จะนำสารสนเทศส่วนนี้ไปประมวลร่วมกับความเป็นไปได้ของพารามิเตอร์ในการแจกแจงความน่าจะเป็นเบื้องต้นก่อน

กำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์

การกำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์ \(\theta\) สามารถทำได้หลายลักษณะ ทั้งแบบไม่ต่อเนื่อง และแบบต่อเนื่อง ขึ้นอยู่กับขอบเขตของการวิเคราะห์ ข้อมูลในอดีตที่ใช้สนับสนุน และธรรมชาติของพารามิเตอร์นั้น ในตัวอย่างข้างต้นพารามิเตอร์ \(\theta\) ต้องมีค่าอยู่บนช่วง \([0,1]\)

เพื่อให้ตัวอย่างนี้ทำความเข้าใจได้ง่าย จึงเลือกกำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์ \(\theta\) เป็นการแจกแจงแบบไม่ต่อเนื่อง โดยมีค่าที่เป็นไปได้เท่ากับ 0.0, 0.1, 0.2, 0.3, …., 1.0 และมีการแจกแจงของความน่าจะเป็นก่อนหน้าดังรูป

# prior distribution
theta<-seq(0,1,0.1)
p.theta<-c(0.01,0.04,0.075,0.1,0.15,0.25,0.15,0.1,0.075,0.04,0.01)
plot(theta,p.theta, 
     xlab=expression(theta), 
    ylab=expression(p(theta)), 
    type="h", lwd=2, col="orange",
    main="Prior Distribution")
points(theta,p.theta, type="p",pch=16, cex=2, col="orange")

การประมาณการแจกแจงความน่าจะเป็นภายหลัง

จากการทดลองโยนเหรียญจำนวน 10 ครั้งพบว่า ออกหน้าหัวจำนวน 6 ครั้ง ดังนั้นจะได้ว่า

\(p(\bf{y}|\theta)=\theta^{6}(1-\theta)^4\)

จากขอบเขตของ \(theta\) ที่กำหนดในการแจกแจงความน่าจะเป็นเบื้องต้นจะได้ว่าฟังก์ชันภาวะความควรจะเป็น มีลักษณะดังรูป

# likelihood function
L<-theta^6*(1-theta)^4
plot(theta,L, 
     xlab=expression(theta), 
    ylab="p(y|theta)",
    type="h", lwd=2, col="orange",
    main="Likelihood Function")
points(theta,L, type="p",pch=16, cex=2, col="orange")

จากทฤษฎีบทของเบส์เราสามารถรวมสารสนเทศจากการแจกแจงความน่าจะเป็นเบื้องต้น กับฟังก์ชันภาวะความควรจะเป็นได้ดังนี้

# marginal y
p.y<-sum(p.theta*L)

par(mfrow=c(1,2))
# posterior of theta given y
post.theta<-L*p.theta/p.y
plot(theta,post.theta, 
     xlab=expression(theta), 
    ylab="p(theta|y)",
    type="h", lwd=2, col="blue",
    ylim=c(0,0.5),
    main="Posterior Distribution")
points(theta,post.theta, type="p",pch=16, cex=2, col="blue")

# prior and posterior in the same plot
plot(theta,p.theta, 
     xlab=expression(theta), 
    ylab="Probability", 
    type="l", lwd=2, col="orange",
    ylim=c(0,0.5))
points(theta,p.theta, type="p",pch=16, cex=2, col="orange")

points(theta,post.theta, 
    type="l", lwd=2, col="blue")
points(theta,post.theta, type="p",pch=16, cex=2, col="blue")
legend(-0.05,0.48, legend=c("Prior","Posterior"), col=c("orange","blue"), lty=1, bty="n", lwd=2)

การอนุมานเชิงสถิติจากการแจกแจงความน่าจะเป็นภายหลัง

เมื่อได้การแจกแจงความน่าจะเป็นภายหลังที่ต้องการแล้ว ผู้วิเคราะห์สามารถวิเคราะห์การแจกแจงดังกล่าวเพื่ออนุมานเกี่ยวกับพารามิเตอร์ความลำเอียงที่ต้องการศึกษาได้ โดยสามารถทำได้หลายลักษณะ ทั้งการประมาณค่าแบบจุด การประมาณค่าแบบช่วง และการทดสอบสมมุติฐาน รายละเอียดดังนี้

(1) การประมาณค่าแบบจุด

การประมาณค่าแบบจุดสามารถทำได้โดยใช้สถิติพื้นฐานสรุปแนวโน้มสู่ส่วนกลางของค่าพารามิเตอร์ จากการแจกแจงความน่าจะเป็นภายหลัง เช่น ค่าเฉลี่ย มัธยฐาน และฐานนิยม นอกจากนี้ยังสามารถคำนวณค่าส่วนเบี่ยงเบนมาตรฐานเพื่อประเมินระดับความน่าเชื่อถือของค่าประมาณแบบจุดดังกล่าวได้อีกด้วย จากการแจกแจงความน่าจะเป็นภายหลังข้างต้น จะได้ว่า

# mean = expected value of theta
mean.theta<-sum(post.theta*theta)
sd.theta<-sqrt(sum(post.theta*(theta-mean.theta)^2))
mean.theta

## [1] 0.5540441

sd.theta

## [1] 0.1145753

(2) การประมาณค่าแบบช่วง

ช่วงการประมาณที่ได้จากวิธีแบบเบส์จะเรียกว่า ช่วงความน่าเชื่อถือ (credible interval) การประมาณค่าพารามิเตอร์แบบช่วงจากการแจกแจงความน่าจะเป็นภายหลังสามารถทำได้หลายวิธีการ วิธีการง่าย ๆ คือการใช้ค่าเฉลี่ยบวกลบด้วยส่วนเบี่ยงเบนมาตรฐานของพารามิเตอร์ อีกวิธีการหนึ่งที่มีประสิทธิภาพมากกว่าคือการหาช่วงแบบ highest density interval (HDI) กล่าวคือเป็นช่วงการประมาณที่มีความหนาแน่นมากที่สุดที่ทำให้ค่าความน่าจะเป็นของพารามิเตอร์ภายในช่วงดังกล่าวเท่ากับค่าที่กำหนดเช่น ช่วง 95% HDI ของพารามิเตอร์ \(\theta\) ในข้างต้นจะมีค่าประมาณ \([0.4,0.8]\)

ช่วงดังกล่าวมีความแตกต่างจาก 95% confidence interval อย่างไร??

(3) การทดสอบสมมุติฐาน

การทดสอบสมมุติฐานด้วยวิธีการแบบเบส์มีความยืดหยุ่นสูงมาก และสามารถทำได้หลายวิธีการ วิธีการแรกเรียกว่า maximum a posteriori (MAP) test มีรายละเอียดดังนี้

กำหนดให้ \(H_0\) และ \(H_1\) เป็นสมมุติฐานที่ต้องการเปรียบเทียบ ในกรณีนี้จะเห็นว่ามีความไม่แน่นอนว่าสมมุติฐานใดเป็นสมมุติฐานที่ถูกต้องกันแน่ ดังนั้นจึงสามารถกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าของสมมุติฐานทั้งสองได้ ดังนี้

\(p(H_0) = p_0\) และ \(p(H_1)=p_1\) โดยที่ \(p_0+p_1=1\)

และกำหนดให้ \(p(\bf{y}|H_0)\) กับ \(p(\bf{y}|H_1)\) คือฟังก์ชันภาวะความควรจะเป็นบนสมมุติฐานทั้งสอง ซึ่งจากทฤษฎีบทของเบส์จะได้ว่า

\(p(H_0|\bf{y})=\frac{p(\bf{y}|H_0)p(H_0)}{p(\bf{y})}\)

\(p(H_1|\bf{y})=\frac{p(\bf{y}|H_1)p(H_1)}{p(\bf{y})}\)

การตัดสินใจว่าควรเลือกเชื่อสมมุติฐานใด สามารถทำได้โดยดูจากอัตราส่วนที่เรียกว่า posterior odd ซึ่งมีค่าเท่ากับ

\(\frac{p(H_0|\bf{y})}{p(H_1|\bf{y})}=\frac{p(\bf{y}|H_0)p(H_0)}{p(\bf{y}|H_1)p(H_1)}\)

โดยหากอัตราส่วนข้างต้นมีค่ามากกว่า 1 แสดงว่า สมมุติฐาน \(H_0\) มีแนวโน้มน่าเชื่อถือมากกว่า ในทางกลับกันหากอัตราส่วนดังกล่าวมีค่าน้อยกว่า 1 แสดงว่าสมมุติฐาน \(H_1\) มีแนวโน้มน่าเชื่อถือมากกว่า

จากตัวอย่างโยนเหรียญ หากต้องการทดสอบว่าเหรียญมีความเที่ยงตรงหรือไม่ อาจกำหนดสมมุติฐานเป็นดังนี้

\(H_0: \theta = 0.5\) vs \(H_1: \theta>0.5\)

prior.H<-c(0.5,0.5) # prior for H0 and H1 respectively.
# calculate likelihood
L.H0<-0.5^6*(1-0.5)^4
theta.H1<-theta[theta>0.5]
L.H1<-sum(theta.H1^6*(1-theta.H1)^4)

# calculate posterior odd
(prior.H[1]*L.H0)/(prior.H[2]*L.H1)

## [1] 0.3727444

อีกวิธีการหนึ่งที่สามารถทำได้เรียกว่า minimum cost hypothesis test รายละเอียดมีดังนี้

กำหนดให้

• \(C_{10}\) คือ cost ของการเลือก \(H1\) โดยที่ \(H_0\) ถูกต้อง (cost ของการเกิด type I error)
• \(C_{01}\) คือ cost ของการเลือก \(H_0\) โดยที่ \(H_1\) ถูกต้อง (cost ของการเกิด type II error)

เกณฑ์การพิจารณาจะใช้ค่าความเสี่ยงภายหลัง (posterior risk) ดังนี้

\(\frac{p(H_0|\bf{y})C_{10}}{p(H_1|\bf{y})C_{01}}\)

จากตัวอย่างโยนเหรียญสมมุติว่า \(C_{10}=5C_{01}\)

c10<-5
c01<-1

# posterior risk of accepting H0
(prior.H[2]*L.H1)*c01

## [1] 0.001309962

# posterior risk of accepting H1
(prior.H[1]*L.H0)*c10

## [1] 0.002441406

# risk ratio
(prior.H[1]*L.H0)*c10/(prior.H[2]*L.H1)*c01

## [1] 1.863722

อีกวิธีการที่ง่ายกว่าสองวิธีการแรกมากคือ วิธี ROPE ซึ่งเป็นการเปรียบเทียบกันระหว่างช่วงความน่าเชื่อถือแบบ HDI กับช่วง ROPE ที่เรียกว่า region of practical equivalence ซึ่งเป็นช่วงที่ยอมรับได้ในทางปฏิบัติ (โดยผู้วิเคราะห์) ว่าหากพารามิเตอร์ยังมีค่าอยู่ในช่วง ROPE ดังกล่าวอยู่ แสดงว่ายังไม่แตกต่างจาก \(H_0\) โดยหากช่วง HDI ครอบคลุมช่วง ROPE เราจะสามารถยอมรับ \(H_0\) ได้ แต่หากช่วง HDI ไม่มีส่วนใดที่ครอบคลุมช่วง ROPE เลย แสดงว่าไม่ช่วงเชื่อถือ \(H_0\) อย่างยิ่ง

การคัดเลือกโมเดล (model selection)

การคัดเลือกโมเดลเป็นวิธีการหนึ่งที่สามารถใช้ในการทดสอบมสมมุติฐานของผู้วิเคราะห์ รวมทั้งใช้เปรียบเทียบและคัดเลือกโมเดลที่เหมาะสมที่สุดสำหรับอธิบายข้อมูล สถิติตัวหนึ่งที่มีความสำคัญคือ bayes factor วัตถุประสงค์ของ bayes factor ถูกใช้เพื่อวัดว่าข้อมูลตัวอย่างที่เก็บรวบรวมมมานี้สนับสนุนโมเดลหรือสมมุติฐานตัวใดมากกว่ากัน bayes factor มีค่าเท่ากับอัตราส่วนระหว่างค่าภาวะความควรจะเป็นของสมมุติฐานสองตัวบนชุดข้อมูลตัวอย่างเดียวกัน สมมุติให้ \(H_1\) คือโมเดลตามสมมุติฐานที่ 1 และ \(H_2\) คือโมเดลตามสมมุติฐานที่ 2 และ \(\bf{x}\) คือชุดข้อมูลตัวอย่าง จะได้ว่าสถิติ bayes factor (BF) มีค่าเท่ากับ

\(BF = \frac{p(\bf{x}|H_1)}{p(\bf{x}|H_2)}=\frac{p(H_1|\bf{x})}{p(H_2|\bf{x})}\times \frac{p(H_2)}{p(H_1)}\)

จากนิยามข้างต้นจะเห็นว่า bayes factor ตามนิยามดังกล่าวสามารถเขียนใหม่ในรูปผลคูณระหว่างอัตราส่วนสองตัว ได้แก่ อัตราส่วนระหว่างการแจกแจงความน่าจะเป็นภายหลังกับอัตราส่วนระหว่างการแจกแจงความน่าจะเป็นก่อนหน้าของสมมุติฐานที่ 1 ต่อสมมุติฐานที่ 2

จากนิยามนี้สามารถแปลผลได้ว่าหาก \(BF\) มีค่ามากกว่า 1 นั่นหมายถึงข้อมูลที่มีสนับสนุนโมเดลตามสมมุติฐาน \(H_1\) มากกว่า ในทางกลับกันหาก \(BF\) มีค่าน้อยกว่า 1 นั่นหมายความว่าข้อมูลที่มีสนับสนุนโมเดลตามสมมุติฐาน \(H_2\) มากกว่า Harold Jeffery (1998) ได้เสนอสเกลการแปลความหมายสถิติ bayes factor ไว้ดังนี้

การคำนวณและการใช้งาน bayes factor จะกล่าวถึงในสัปดาห์ที่สองของหัวข้อนี้

นอกจากสถิติ bayes factor แล้วยังมีสถิติอีกหลายตัวที่สามารถใช้เปรียบเทียบและคัดเลือกโมเดลแบบเบส์ได้ เช่น deviance information criterion (DIC), widely applicable information criterion (WAIC) และ leave-one-out cross-validation (LOO) สถิติทั้งหมดที่กล่าวมานี้ใช้สะท้อนความแตกต่างกันระหว่างข้อมูลตัวอย่าง \(y\) กับค่าทำนายที่ generate จากโมเดลการวิเคราะห์ \(\hat{y}\) ค่าที่ได้จากสถิติดังกล่าวจึงใช้สะท้อนความคลาดเคลื่อนของแต่ละโมเดลเมื่อเปรียบเทียบกับข้อมูลจริง โมเดลที่มีค่าสถิติดังกล่าวต่ำกว่าจะเป็นโมเดลที่มีความเหมาะสมมากกว่า

การแจกแจงความน่าจะเป็นก่อนหน้า

การแจกแจงความน่าจะเป็นก่อนหน้าใช้สะท้อนความเชื่อเบื้องต้นในพารามิเตอร์ต่าง ๆ ของโมเดล ความเชื่อดังกล่าวมาได้จากทั้งผลการศึกษาในอดีต หลักเหตุผล และประสบการณ์ของนักวิจัย นอกจากนี้ยังควรต้องพิจารณาธรรมชาติหรือค่าที่เป็นไปได้ของพารามิเตอร์ประกอบด้วย

โดยทั่วไปการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าสามารถทำได้ 2 ลักษณะ ลักษณะแรกเรียกว่า การแจกแจงความน่าจะเป็นก่อนหน้าแบบไม่ให้สารสนเทศ (non-informative prior distribution) และการแจกแจงความน่าจะเป็นก่อนหน้าแบบให้สารสนเทศ (informative prior distribution)

การแจกแจงความน่าจะเป็นก่อนหน้าของพารามิเตอร์อิทธิพลคงที่สามารถเลือกใช้ได้หลายการแจกแจง เช่น การแจกแจงแบบ uniform ที่เป็น noninformative prior การแจกแจงแบบปกติ ที่เป็นไปได้ทั้ง noninformative และ informative prior ขึ้นอยู่กับการกำหนดพารามิเตอร์ค่าเฉลี่ย (mean) และค่าความถูกต้อง (precision) โดยที่ค่าความถูกต้องดังกล่าวเป็นส่วนกลับของพารามิเตอร์ความแปรปรวน (variance) ดังนี้ precision = 1/variance

จากตัวอย่างข้างต้นกำหนดให้

\(\beta_0 \sim N(0,10^{-4})\)

\(\beta_1 \sim N(0,10^{-4})\)

การกำหนดข้างต้นคือการแจกแจงความน่าจะเป็นแบบปกติที่มีค่าเฉลี่ยเท่ากับ 0 และความแปรปรวนเท่ากับ 10000 ซึ่งเรียกได้ว่าเป็นการแจกแจงความน่าจะเป็นก่อนหน้าแบบไม่ให้สารสนเทศ รูปด้่านล่างแสดงลักษณะการแจกแจงแบบปกติข้างต้น จากรูปจะเห็นว่าการแจกแจงดังกล่าวมีขอบเขตที่กว้างพารามิเตอร์ส่วนใหญ่จะมีค่าตกอยู่ในช่วง -200 ถึง 200

การแจกแจงความน่าจะเป็นก่อนหน้าของส่วนเบี่ยงเบนมาตรฐานหรือความแปรปรวน มีความจำเป็นต้องเลือกการแจกแจงที่ domain ของการแจกแจงมีค่าไม่ติดลบ โดยปกติมักกำหนดให้ความแปรปรวนดังกล่าวมีการแจกแจงแบบแกมมาผกผัน (inverse gamma distribution)

จากตัวอย่างข้างต้น กำหนดให้

\(sigma^2 \sim IG(a,b)\)

โดยที่ \(IG\) คือการแจกแจงแบบ inverse-gamma

Note: ในเชิงทฤษฎีหากกำหนดให้ \(sigma^2 \sim IG(a,b)\) จะได้ว่า

\(p(\sigma^2|a,b) \propto (\sigma^2)^{-a-1}exp(\frac{-b}{\sigma^2})\)

หากกำหนดให้ \(a \rightarrow 0\) และ \(b \rightarrow 0\) แล้วฟังก์ชันความหนาแน่นในข้างต้นจะเขียนใหม่ได้เป็น

\(p(\sigma^2|a,b) \propto \frac{1}{\sigma^2}\) เรียกว่า Jeffery prior ซึ่งเป็น noninformative prior ที่ใช้ได้ตัวหนึ่งของพารามิเตอร์ความแปรปรวน

ดังนั้นในกรณีนี้จึงกำหนดให้

\(sigma^2 \sim IG(0.0001,0.0001)\)

รูปด้านล่างแสดงรายการของการแจกแจงความน่าจะเป็นที่สามารถกำหนดได้ใน BUGS และ JAGS

จากการระบุโมเดลข้างต้นทำให้สามารถเขียน syntax เพื่อระบุโมเดลใน JAGS ได้ดังนี้

model{
  
  for(i in 1:n)
  {
    # likelihood function
    y[i]~dnorm(mu[i],tau) #where tau is precision = 1/sigma2
    
    # deteministic model
    mu[i]<-b0+b1x[i] 
  }
  
  # prior distribution
  b0~dnorm(0,10^-4)
  b1~dnorm(0,10^-4)
  
  tau~gamma(0.0001,0.0001)
  sigma2<-1/tau # deteministic model
}

การประมวลผลด้วย JAGS

การประมวลผลด้วย JAGS มี 4 ขั้นตอนหลักได้แก่

(1) เตรียมข้อมูล

(2) ระบุโมเดล

(3) ส่งข้อมูลและโมเดลไปประมวลผลบน JAGS

(4) ตรวจสอบการลู่เข้าของ MCMC

(5) วิเคราะห์และแปลผล

modelString<-"model{
  
  for(i in 1:n)
  {
    # likelihood function
    y[i]~dnorm(mu[i],tau) #where tau is precision = 1/sigma2
    
    # deteministic model
    mu[i]<-b0+b1*x[i] 
  }
  
  # prior distribution
  b0~dnorm(0,10^-4)
  b1~dnorm(0,10^-4)
  
  tau~dgamma(0.0001,0.0001)
  sigma2<-1/tau # deteministic model
}" #Wrap model syntax into String object

writeLines(modelString, con="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/simReg.txt") # write modelString to model.txt

ATT<-rnorm(100,30,10)
ACH<-30+0.55*ATT+rnorm(100,0,5)
dat.reg<-data.frame(ATT,ACH)

head(dat.reg)

##         ATT      ACH
## 1 32.133899 52.34851
## 2  8.208871 37.08404
## 3 37.542008 54.13731
## 4 13.544627 38.93896
## 5 35.487326 44.67005
## 6 42.563916 50.66076

par(mar=c(4,4,1,1))
plot(ATT,ACH, pch=16, xlab="ATT", ylab="ACH")

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

library(rjags)
dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
initsList<-list(b0=1, b1=1, tau=1)

model<-jags.model(file = "simReg.txt",
                  data = dataList,
                  inits = initsList,
                  n.chains = 3)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 3
##    Total graph size: 412
## 
## Initializing model

update(model, n.iter=2000) # n.burnin
samples<-coda.samples(model,
                      variable.names=c("b0","b1","sigma2"),
                      n.iter = 10000,
                      thin=1)
head(samples)

## [[1]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 2001 
## End = 2007 
## Thinning interval = 1 
##            b0        b1   sigma2
## [1,] 30.78025 0.5202950 21.90599
## [2,] 30.60087 0.4840719 24.14094
## [3,] 31.66543 0.4708754 37.18908
## [4,] 31.32383 0.4505433 21.38748
## [5,] 31.89650 0.4600205 24.17811
## [6,] 31.78018 0.4470285 30.98406
## [7,] 32.44455 0.4417329 24.82733
## 
## [[2]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 2001 
## End = 2007 
## Thinning interval = 1 
##            b0        b1   sigma2
## [1,] 29.98931 0.4762332 30.62127
## [2,] 29.89494 0.5395243 21.85398
## [3,] 29.32648 0.5483189 30.41757
## [4,] 29.36315 0.5429597 31.03104
## [5,] 29.14981 0.5379856 22.40325
## [6,] 30.14993 0.5327550 26.61984
## [7,] 29.80438 0.5197910 23.47122
## 
## [[3]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 2001 
## End = 2007 
## Thinning interval = 1 
##            b0        b1   sigma2
## [1,] 31.81481 0.4581798 24.36047
## [2,] 31.77167 0.4908677 20.95915
## [3,] 31.31646 0.4657395 25.43489
## [4,] 31.52530 0.4829907 30.82938
## [5,] 31.05918 0.4862660 22.68679
## [6,] 30.97493 0.5190658 24.36354
## [7,] 30.30192 0.5135753 24.93864
## 
## attr(,"class")
## [1] "mcmc.list"

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

library(runjags)
dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
initsList<-list(b0=1, b1=1, tau=1)

runjag.output<-run.jags(method="parallel",
                model="simReg.txt",
                monitor=c("b0","b1","sigma2"),
                data=dataList,
                inits=initsList,
                n.chains=3,
                burnin=2000,
                sample=10000, #n.iter
                thin=1,
                summarise=FALSE,
                plots=FALSE)
  
samples<-as.mcmc.list(runjag.output)
head(samples)
plot(samples)

library(MCMCvis)
MCMCtrace(samples, params=c("b0","b1","sigma2"), pdf=F)

MCMCsummary(samples)

##              mean         sd       2.5%        50%      97.5% Rhat n.eff
## b0     30.0629002 1.61239595 26.9107602 30.0698484 33.2015106    1  1633
## b1      0.5213053 0.05160855  0.4217729  0.5211085  0.6222152    1  1607
## sigma2 26.9798622 3.93882192 20.3574546 26.6128567 35.6623178    1 27719

summary(samples)

## 
## Iterations = 2001:12000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD Naive SE Time-series SE
## b0     30.0629 1.61240 0.009309       0.039919
## b1      0.5213 0.05161 0.000298       0.001288
## sigma2 26.9799 3.93882 0.022741       0.023665
## 
## 2. Quantiles for each variable:
## 
##           2.5%     25%     50%    75%   97.5%
## b0     26.9108 28.9716 30.0698 31.166 33.2015
## b1      0.4218  0.4863  0.5211  0.556  0.6222
## sigma2 20.3575 24.1980 26.6129 29.370 35.6623

autocorr.diag(samples)

##                b0         b1      sigma2
## Lag 0  1.00000000 1.00000000 1.000000000
## Lag 1  0.89674469 0.89729301 0.018382972
## Lag 5  0.58881031 0.58970557 0.002997731
## Lag 10 0.34899078 0.35034313 0.004552666
## Lag 50 0.03019034 0.03073331 0.004531337

autocorr.plot(samples[[1]])

การอนุมานเชิงสถิติแบบเบส์

จากผลการวินิจฉัย MCMC ในข้างต้น เราทำการประมวลผลใหม่ดังนี้

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

#dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
#initsList<-list(b0=1, b1=1, tau=1)

model<-jags.model(file = "simReg.txt",
                  data = dataList,
                  inits = initsList,
                  n.chains = 3)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 3
##    Total graph size: 412
## 
## Initializing model

update(model, n.iter=2000) # n.burnin
samples<-coda.samples(model,
                      variable.names=c("b0","b1","sigma2"),
                      n.iter = 30000,
                      thin=15)
MCMCsummary(samples)

##              mean         sd       2.5%        50%      97.5% Rhat n.eff
## b0     29.9358012 1.61122865 26.7793446 29.9515443 33.0316408    1  4086
## b1      0.5251646 0.05136625  0.4260559  0.5245155  0.6249333    1  4067
## sigma2 26.9815883 3.93464366 20.3197188 26.6338984 35.8743312    1  6000

MCMCtrace(samples, params=c("b0","b1","sigma2"), pdf=F)

autocorr.diag(samples)

##                   b0           b1       sigma2
## Lag 0    1.000000000 1.0000000000  1.000000000
## Lag 15   0.189732427 0.1917814786 -0.003137443
## Lag 75   0.026335394 0.0191889064  0.007447446
## Lag 150  0.006470748 0.0103972916 -0.007191825
## Lag 750 -0.003792101 0.0003762561  0.025105044

autocorr.plot(samples[[1]])

summary(samples)

## 
## Iterations = 2015:32000
## Thinning interval = 15 
## Number of chains = 3 
## Sample size per chain = 2000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD  Naive SE Time-series SE
## b0     29.9358 1.61123 0.0208009      0.0252067
## b1      0.5252 0.05137 0.0006631      0.0008054
## sigma2 26.9816 3.93464 0.0507960      0.0507991
## 
## 2. Quantiles for each variable:
## 
##           2.5%   25%     50%    75%   97.5%
## b0     26.7793 28.84 29.9515 31.033 33.0316
## b1      0.4261  0.49  0.5245  0.561  0.6249
## sigma2 20.3197 24.21 26.6339 29.332 35.8743

#library(googledrive)
setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel")
#drive_download("https://drive.google.com/file/d/0B3lh4V2Mrl14SF94US1YTkNVaHM/view?resourcekey=0-bYI7GCSuypi3a5aFt0O-qg")
source("DBDA2E-utilities.R") # (Kruschke, 2015)

## 
## *********************************************************************
## Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:
## A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
## *********************************************************************

plotPost(samples[,"b1"], compVal=0, cenTend="median", ROPE=c(0,0.1), credMass=0.95)

##                  ESS      mean    median      mode hdiMass    hdiLow   hdiHigh
## Param. Val. 4069.857 0.5251646 0.5245155 0.5166871    0.95 0.4260469 0.6249327
##             compVal pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## Param. Val.       0          1       0      0.1       0       0       1

plotPost(samples[,"b1"], compVal=0.6, cenTend="median", ROPE=c(0.55,0.65), credMass=0.95)

##                  ESS      mean    median      mode hdiMass    hdiLow   hdiHigh
## Param. Val. 4069.857 0.5251646 0.5245155 0.5166871    0.95 0.4260469 0.6249327
##             compVal pGtCompVal ROPElow ROPEhigh   pLtROPE   pInROPE     pGtROPE
## Param. Val.     0.6 0.07133333    0.55     0.65 0.6831667 0.3091667 0.007666667

ตัวอย่างโยนเหรียญ (รอบที่ 2)

• จากตัวอย่างที่ 1 จะได้ว่าโมเดลของค่าสังเกตเป็นฟังก์ชันความน่าจะเป็นแบบ Bernuolli ที่มีพารามิเตอร์เป็น \(\theta\) เขียนแทนด้วย

• อย่างไรก็ตามในกรณีนี้ผู้วิเคราะห์ต้องการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าแบบต่อเนื่องให้กับพารามิเตอร์ความลำเอียง \(\theta\) เนื่องจากพารามิเตอร์ดังกล่าวต้องมีค่าอยู่บนช่วง \([0,1]\) การแจกแจงความน่าจะเป็นที่เหมาะสมจึงเป็นการแจกแจงแบบบีต้า (beta distribution) ที่มีพารามิเตอร์ตำแหน่ง \(a\) และพารามิเตอร์สเกล \(b\) ดังนั้น

พารามิเตอร์ทั้งสองตัวของการแจกแจงแบบบีต้าทำให้รูปทรงการแจกแจงมีความแตกต่างกัน ซึ่งสามารถใช้สะท้อนความเชื่อเกี่ยวกับความลำเอียงของเหรียญที่แตกต่างกันได้เป็นอย่างดี รูปด้านล่างแสดงตัวอย่างโค้งความหนาแน่นของการแจกแจงแบบบีต้า เมื่อกำหนดพารามิเตอร์ต่าง ๆ

Kruschke (2012)

จากการกำหนดในข้างต้นจะสามารถเขียน syntax ของโมเดลในโปรแรกม JAGS ได้ดังนี้

model{

# likelihood part    
for (i in 1:n)
{
  y[i]~dbern(theta)
  
}
  
# prior part
  theta~beta(a,b)

# deterministic part
  a<-3
  b<-2
}

การประมวลผลโมเดลข้างต้นด้วยโปรแกรม JAGS บน R จำเป็นต้องมีการ wrap up คำสั่งข้างต้นให้เป็น text file เพื่อส่งออกไปยังโปรแกรม R โดยเขียนคำสั่งดังนี้

coin_model="

model{

# likelihood part    
for (i in 1:n)
{
  y[i]~dbern(theta)
  
}
  
# prior part
  theta~dbeta(a,b)

# deterministic part
  a<-3
  b<-2
}

"
writeLines(coin_model, con="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/coin_model.txt")

ขั้นตอนถัดมาคือการประมาณการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ \(\theta\) โดยใช้เทคนิค MCMC ผ่านโปรแกรม JAGS ซึ่งสามารถสั่งงานได้ผ่านฟังก์ชัน jag.model() โดยฟังก์ชันนี้มีอาร์กิวเมนท์ที่สำคัญที่เกี่ยวข้องกับการระบุคุณลักษณะของลูกโซ่มาร์คอฟที่ต้องการสร้างขึ้น ได้แก่

• ข้อมูลค่าสังเกต ในที่นี้สมมุติว่าทำการทดลอง 20 ครั้งได้ผลเป็นดังนี้

y<-rbinom(20,1,0.65)
data.list<-list(y=y, n=20)

• n.chains

library(rjags)
fit<-jags.model(file="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/coin_model.txt", data=data.list,
                n.chains = 3)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 20
##    Unobserved stochastic nodes: 1
##    Total graph size: 24
## 
## Initializing model

samples<-coda.samples(fit,variable.names=c("theta"),n.iter=5000)

plot(samples)

autocorr.diag(samples)

##              theta
## Lag 0   1.00000000
## Lag 1   0.01582687
## Lag 5  -0.01269683
## Lag 10  0.01048772
## Lag 50 -0.01040355

summary(samples)

## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##      0.6796919      0.0908805      0.0007420      0.0007499 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.4931 0.6189 0.6847 0.7465 0.8412

เหรียญเที่ยงตรงหรือไม่?? —> กำหนด ROPE = [0.48,0.52]

plotPost(samples[,"theta"], compVal=0.5, cenTend="median", ROPE=c(0.48,0.1), credMass=0.95)

##                  ESS      mean    median      mode hdiMass    hdiLow   hdiHigh
## Param. Val. 14532.27 0.6796919 0.6847037 0.7075729    0.95 0.5029832 0.8492614
##             compVal pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## Param. Val.     0.5     0.9702    0.48      0.1  0.0186       0       1

One-Way ANOVA with random effect

null model หรือโมเดลการวิเคราะห์ความแปรปรวนแบบอิทธิพลสุ่มมีสมการดังนี้

level-1 model: \(y_{ij}=\beta_{0j}+\epsilon_{ij}\)

level-2 model: \(\beta_{0j}=\gamma_{00}+u_{0j}\)

combined model: \(y_{ij}=\gamma_{00}+u_{0j}+\epsilon_{ij}\)

โดยที่ตัวแปรตามคือ mathach ในชุดข้อมูล hsb หากดำเนินการวิเคราะห์ข้อมูลด้วยวิธีการแบบดั้งเดิม จะได้ผลการวิเคราะห์เป็นดังนี้

library(foreign)
dat1<-read.spss("/Users/siwachoat/Downloads/hsb1.sav", to.data.frame = TRUE)
dat2<-read.spss("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/เอกสารประกอบการสอน/2756714/HLM/hsb2.sav", to.data.frame=TRUE)
dat<-merge(dat1,dat2, by="schoolid")
library(lme4)

## Loading required package: Matrix

library(lmerTest)

## 
## Attaching package: 'lmerTest'

## The following object is masked from 'package:lme4':
## 
##     lmer

## The following object is masked from 'package:stats':
## 
##     step

fit<-lmer(mathach~1+(1|schoolid), data=dat1)
ranova(fit) #reduce random-effect term

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## mathach ~ (1 | schoolid)
##                npar logLik   AIC    LRT Df Pr(>Chisq)    
## <none>            3 -23558 47123                         
## (1 | schoolid)    2 -24052 48107 986.12  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mathach ~ 1 + (1 | schoolid)
##    Data: dat1
## 
## REML criterion at convergence: 47116.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0631 -0.7539  0.0267  0.7606  2.7426 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schoolid (Intercept)  8.614   2.935   
##  Residual             39.148   6.257   
## Number of obs: 7185, groups:  schoolid, 160
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  12.6370     0.2444 156.6473   51.71   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

performance::icc(fit)

## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.180
##   Conditional ICC: 0.180

# reliability of the sample mean (Raudenbush & Bryk (2002) p.46)
nj<-table(dat1$schoolid)
reliaj<-8.614/(8.614+39.148/nj)
hist(reliaj)

mean(reliaj)

## [1] 0.9013778

#calculate level-1 residual
res1<-dat1$mathach-predict(fit)
res1<-data.frame(schoolid=dat1$schoolid, res1)
#calculate level-2 residual
res2<-coef(fit)$schoolid-12.6370
names(res2)<-"res2"
res2<-data.frame(schoolid=as.numeric(names(table(dat1$schoolid))), res2)


#normality
ggplot(res1)+geom_histogram(aes(x=res1))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(res2)+geom_histogram(aes(x=res2))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

library(ggpubr)
ggqqplot(res1$res1)

ggqqplot(res2$res2)

#homogeneity of variance
res1%>%group_by(schoolid)%>%summarise(var.e=var(res1))%>%summary()

##     schoolid        var.e      
##  Min.   :1224   Min.   :12.54  
##  1st Qu.:3018   1st Qu.:31.27  
##  Median :5512   Median :39.94  
##  Mean   :5310   Mean   :39.75  
##  3rd Qu.:7432   3rd Qu.:47.10  
##  Max.   :9586   Max.   :71.93

res1%>%group_by(schoolid)%>%summarise(var.e=var(res1))%>%ggplot(aes(x=var.e))+geom_density()

หากดำเนินการวิเคราะห์ด้วยวิธีการแบบ bayes โดยใช้โปรแกรม JAGS อาจดำเนินการได้ดังนี้

library(rjags)

#generate model file
null.model<-"model{

for (i in 1:n)
{
y[i]~dnorm(mu[i],inv.sigma2)
mu[i]<-beta0j[schoolid[i]]
}

for (j in tab)
{
beta0j[j]~dnorm(gamma00, inv.tau00)
}

inv.sigma2~dgamma(0.0001,0.0001)
inv.tau00~dgamma(0.0001,0.0001)
gamma00~dnorm(0,10^-4) T(0,) #half-normal distribution

sigma2<-inv.sigma2^(-1)
tau00<-inv.tau00^(-1)

icc<-tau00/(tau00+sigma2)

}"

writeLines(null.model,con="null_model.txt")

datlist<-list(y=dat1$mathach,schoolid=dat1$schoolid, n=length(dat1$mathach),
              tab=as.numeric(names(table(dat1$schoolid))))
fit.null<-jags.model(file="null_model.txt",
                     data=datlist,
                     n.chains=3)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7185
##    Unobserved stochastic nodes: 163
##    Total graph size: 14706
## 
## Initializing model

update(fit.null,n.iter=2000) #n.burnin
output<-coda.samples(fit.null,
                     variable.names=c("gamma00","tau00","sigma2","icc"),
                     n.iter=30000,
                     thin=2)

library(MCMCvis)
MCMCtrace(output, params=c("gamma00","icc"), pdf=F)

MCMCsummary(output)

##               mean         sd       2.5%        50%      97.5% Rhat n.eff
## gamma00 12.6367876 0.24487553 12.1583073 12.6381989 13.1172507    1 45000
## icc      0.1817391 0.01885253  0.1474386  0.1808955  0.2212261    1 40991
## sigma2  39.1627730 0.65410399 37.9016391 39.1535162 40.4598172    1 44151
## tau00    8.7211157 1.10051603  6.7886400  8.6454339 11.0962064    1 40938

summary(output)

## 
## Iterations = 2002:32000
## Thinning interval = 2 
## Number of chains = 3 
## Sample size per chain = 15000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean      SD  Naive SE Time-series SE
## gamma00 12.6368 0.24488 1.154e-03      1.154e-03
## icc      0.1817 0.01885 8.887e-05      9.312e-05
## sigma2  39.1628 0.65410 3.083e-03      3.114e-03
## tau00    8.7211 1.10052 5.188e-03      5.440e-03
## 
## 2. Quantiles for each variable:
## 
##            2.5%     25%     50%     75%   97.5%
## gamma00 12.1583 12.4714 12.6382 12.8033 13.1173
## icc      0.1474  0.1686  0.1809  0.1939  0.2212
## sigma2  37.9016 38.7157 39.1535 39.6046 40.4598
## tau00    6.7886  7.9512  8.6454  9.4048 11.0962

plotPost(output[,"icc"], compVal=0.05, cenTend="median", ROPE=c(0.04,0.06), credMass=0.95,
         xlab="ICC")

##          ESS      mean    median      mode hdiMass    hdiLow   hdiHigh compVal
## ICC 41391.37 0.1817391 0.1808955 0.1794032    0.95 0.1457078 0.2191428    0.05
##     pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## ICC          1    0.04     0.06       0       0       1

library(rjags)

#generate model file
null.model<-"model{

for (i in 1:n)
{
y[i]~dnorm(mu[i],inv.sigma2[schoolid[i]])
mu[i]<-beta0j[schoolid[i]]
}

for (j in tab)
{
beta0j[j]~dnorm(gamma00, inv.tau00)
inv.sigma2[j]~dgamma(0.0001,0.0001)
sigma2[j]<-inv.sigma2[j]^(-1)

}

inv.tau00~dgamma(0.0001,0.0001)
gamma00~dnorm(0,10^-4) T(0,) #half-normal distribution
tau00<-inv.tau00^(-1)
}"

writeLines(null.model,con="null_model.txt")

datlist<-list(y=dat1$mathach,schoolid=dat1$schoolid, n=length(dat1$mathach),
              tab=as.numeric(names(table(dat1$schoolid))))
fit.hetero<-jags.model(file="null_model.txt",
                     data=datlist,
                     n.chains=3)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7185
##    Unobserved stochastic nodes: 322
##    Total graph size: 15022
## 
## Initializing model

update(fit.hetero,n.iter=2000) #n.burnin
output<-coda.samples(fit.hetero,
                     variable.names=c("gamma00","tau00","sigma2"),
                     n.iter=30000,
                     thin=2)

library(MCMCvis)
MCMCsummary(output)

##                  mean         sd      2.5%       50%     97.5% Rhat n.eff
## gamma00      12.65549  0.2514842 12.160433 12.655888  13.14829    1 44242
## sigma2[1224] 60.19157 13.0677194 39.872462 58.449081  90.65491    1 43915
## sigma2[1288] 53.16940 16.5107789 30.066191 50.173312  93.34146    1 48202
## sigma2[1296] 29.98479  6.5156131 19.869552 29.097814  45.04621    1 45000
## sigma2[1308] 41.80058 15.0983901 21.594343 38.774190  79.14488    1 43913
## sigma2[1317] 31.12595  6.7245490 20.631624 30.210267  46.78136    1 45139
## sigma2[1358] 36.93542 10.4336943 21.831424 35.186547  62.27490    1 44588
## sigma2[1374] 75.43277 22.2474597 43.540262 71.621535 129.38767    1 44960
## sigma2[1433] 16.04133  4.1573219  9.836660 15.405604  25.92200    1 45000
## sigma2[1436] 21.78928  4.9608527 14.178401 21.087452  33.52247    1 45000
## sigma2[1461] 51.46226 13.7561788 31.219932 49.250434  84.58928    1 45000
## sigma2[1462] 41.23397  8.1201423 28.301439 40.212561  59.87144    1 45751
## sigma2[1477] 52.86930  9.8575519 37.028408 51.740710  75.52237    1 45000
## sigma2[1499] 41.90864  8.5942796 28.323192 40.814185  61.53317    1 45000
## sigma2[1637] 55.70535 17.0675853 31.773920 52.783438  97.24919    1 43574
## sigma2[1906] 44.22319  9.0872565 29.905604 43.024882  65.20067    1 45000
## sigma2[1909] 40.72442 12.0302599 23.501970 38.643688  69.78552    1 45689
## sigma2[1942] 32.56230  9.4782120 19.002048 30.980386  55.38963    1 45000
## sigma2[1946] 51.27276 12.3699894 32.584266 49.428249  80.57575    1 45000
## sigma2[2030] 41.28092  8.9668134 27.366108 40.057587  62.26450    1 44440
## sigma2[2208] 38.82988  7.4463692 26.997708 37.934845  55.96154    1 43710
## sigma2[2277] 32.03703  6.0690991 22.234330 31.296887  46.02290    1 45000
## sigma2[2305] 25.89295  4.6239434 18.419459 25.367988  36.54755    1 44896
## sigma2[2336] 35.78740  7.7935331 23.687094 34.738324  53.81243    1 44321
## sigma2[2458] 35.37772  6.9451626 24.317726 34.526280  51.29790    1 45465
## sigma2[2467] 47.85396  9.8813976 32.218197 46.618612  70.90522    1 43913
## sigma2[2526] 23.69674  4.6673615 16.306127 23.134159  34.57381    1 44077
## sigma2[2626] 41.11736 10.1278780 25.901219 39.619390  65.11114    1 44599
## sigma2[2629] 27.64507  5.4027537 19.065542 26.954774  39.99894    1 44623
## sigma2[2639] 36.19415  8.4455039 23.273651 35.045227  56.18398    1 45000
## sigma2[2651] 51.86463 12.6833077 32.668867 49.937409  82.07508    1 44651
## sigma2[2655] 41.53026  8.5117493 28.071720 40.401609  61.29319    1 44513
## sigma2[2658] 33.26988  7.4519265 21.808922 32.246284  50.66791    1 44536
## sigma2[2755] 38.57854  8.3794385 25.506217 37.439966  58.06386    1 45493
## sigma2[2768] 57.43159 18.0613354 32.075967 54.181315 102.02909    1 44939
## sigma2[2771] 47.96107  9.5784090 32.878142 46.737742  70.15289    1 46234
## sigma2[2818] 31.87086  7.3520616 20.544201 30.856273  49.23508    1 43681
## sigma2[2917] 52.78165 12.1475424 34.119049 51.037604  81.29278    1 45385
## sigma2[2990] 19.72811  4.2542055 13.135323 19.173486  29.57098    1 45215
## sigma2[2995] 55.84869 12.3698831 36.606600 54.212405  84.90843    1 44571
## sigma2[3013] 50.60847 10.2515669 34.335236 49.336139  74.20117    1 44502
## sigma2[3020] 41.70495  8.0386673 28.873111 40.713629  60.06811    1 45000
## sigma2[3039] 32.19949 11.4030505 16.956329 29.984724  60.32980    1 45247
## sigma2[3088] 37.93164  9.2291605 23.977790 36.545053  59.81383    1 44641
## sigma2[3152] 56.88409 11.6967789 38.496180 55.370426  83.84153    1 45000
## sigma2[3332] 44.77404 10.9695837 28.151362 43.143046  70.76379    1 44435
## sigma2[3351] 52.18636 12.7072259 32.833498 50.331278  82.03890    1 44321
## sigma2[3377] 51.64357 11.5700697 33.674183 50.079745  78.78349    1 44328
## sigma2[3427] 13.15345  2.8222248  8.736623 12.785176  19.80100    1 45390
## sigma2[3498] 25.97002  5.2919369 17.564168 25.296574  38.18994    1 44511
## sigma2[3499] 31.45782  7.6953345 19.825063 30.341658  49.57740    1 44924
## sigma2[3533] 54.56205 11.7384429 36.166382 53.048400  81.93294    1 45000
## sigma2[3610] 35.87454  6.6020964 25.225146 35.109942  50.95708    1 44821
## sigma2[3657] 62.20887 12.8843325 41.801134 60.576876  91.98888    1 45000
## sigma2[3688] 47.29698 10.8148645 30.614391 45.722289  72.96350    1 45000
## sigma2[3705] 40.55032  9.0918732 26.440599 39.296417  61.97958    1 44033
## sigma2[3716] 75.64338 17.8456689 48.458948 73.001054 117.48385    1 44613
## sigma2[3838] 27.07179  5.4739627 18.448212 26.343909  39.78723    1 45517
## sigma2[3881] 45.88496 10.7714342 29.443579 44.326782  71.49261    1 44448
## sigma2[3967] 49.32437 10.1505469 33.285648 48.019856  72.86548    1 45000
## sigma2[3992] 32.48634  6.6253652 22.033897 31.639625  47.82663    1 45000
## sigma2[3999] 63.68147 14.0315848 41.771119 61.796220  96.44179    1 44880
## sigma2[4042] 42.88707  7.8937138 30.045636 41.997704  60.91079    1 44197
## sigma2[4173] 51.35420 11.5316829 33.611976 49.730882  78.11309    1 45000
## sigma2[4223] 43.00422  9.5897565 28.082518 41.727276  65.60563    1 43404
## sigma2[4253] 33.19243  6.4509104 22.836672 32.397621  47.98253    1 45000
## sigma2[4292] 39.93381  7.2753959 28.219372 39.128877  56.52703    1 45000
## sigma2[4325] 37.66066  7.7088590 25.559442 36.704958  55.60960    1 45000
## sigma2[4350] 55.17794 14.7413463 33.332433 52.847586  90.26246    1 45000
## sigma2[4383] 59.97786 18.6782108 33.838398 56.594689 105.22696    1 44550
## sigma2[4410] 40.21179  9.4435115 25.800708 38.877802  62.31101    1 45000
## sigma2[4420] 45.92915 12.4932123 27.541042 43.901113  75.81138    1 45889
## sigma2[4458] 20.98883  4.5233121 13.885132 20.393734  31.58704    1 45000
## sigma2[4511] 36.62661  7.1252640 25.238617 35.748666  53.01316    1 45359
## sigma2[4523] 32.11103  7.0526311 21.151461 31.177179  48.55198    1 45000
## sigma2[4530] 32.08005  5.9805556 22.532385 31.370734  45.65906    1 45000
## sigma2[4642] 38.35838  7.2570457 26.698959 37.467132  54.92376    1 45819
## sigma2[4868] 31.26781  8.1631762 19.056097 30.049310  50.91516    1 45353
## sigma2[4931] 29.39309  5.6848249 20.333775 28.712575  42.46148    1 45000
## sigma2[5192] 53.85589 15.8712952 31.197746 51.206495  92.15191    1 46876
## sigma2[5404] 37.61934  7.3585860 25.830609 36.738920  54.61382    1 45374
## sigma2[5619] 54.61691  9.9104725 38.674234 53.422470  77.17139    1 44561
## sigma2[5640] 52.25551 10.2331550 35.910943 50.940000  76.04966    1 45000
## sigma2[5650] 43.82994  9.7287861 28.771032 42.482415  66.36990    1 45000
## sigma2[5667] 44.81248  8.4603252 31.335715 43.782942  63.99972    1 45000
## sigma2[5720] 33.67890  6.8762154 22.873839 32.797800  49.50510    1 45000
## sigma2[5761] 44.46468  9.1236501 30.022948 43.347238  65.35517    1 45000
## sigma2[5762] 26.76025  6.7917153 16.587019 25.716064  43.04878    1 45000
## sigma2[5783] 45.67161 13.1608142 26.716212 43.507649  77.71426    1 44609
## sigma2[5815] 34.47724 10.9922264 19.261471 32.499961  61.56169    1 44300
## sigma2[5819] 50.57602 10.6391394 33.959650 49.133487  75.24668    1 45000
## sigma2[5838] 33.30134  9.1853551 19.842431 31.875482  55.33488    1 44572
## sigma2[5937] 26.60941  7.7254317 15.450716 25.331708  45.16086    1 44483
## sigma2[6074] 41.08376  8.1215975 28.156115 40.061159  59.88008    1 44503
## sigma2[6089] 45.04725 12.0059738 27.268417 43.228287  73.90989    1 44506
## sigma2[6144] 44.61165 10.2848589 28.876943 43.153325  68.94483    1 44958
## sigma2[6170] 76.39065 26.8142096 40.705098 71.070778 142.82118    1 45000
## sigma2[6291] 46.19685 11.9690269 28.443755 44.267533  74.81003    1 44036
## sigma2[6366] 41.06237  7.9854519 28.318932 40.062493  59.29037    1 46221
## sigma2[6397] 49.94028  9.4773093 34.827726 48.838631  71.70610    1 44977
## sigma2[6415] 47.60867  9.6312035 32.341571 46.377283  69.74527    1 44451
## sigma2[6443] 37.55948 10.5376479 22.266250 35.815627  62.88573    1 45000
## sigma2[6464] 31.15446  9.0638023 18.119497 29.583505  53.03124    1 45000
## sigma2[6469] 23.67562  4.6462931 16.264451 23.107028  34.35422    1 45000
## sigma2[6484] 58.25340 14.9735380 35.958298 56.000296  93.57001    1 45000
## sigma2[6578] 41.36688  8.1544071 28.339339 40.377228  60.09285    1 46352
## sigma2[6600] 48.59893  9.5948826 33.339704 47.421753  70.60684    1 45000
## sigma2[6808] 51.30961 11.5608818 33.441305 49.707366  78.29894    1 45385
## sigma2[6816] 22.01430  4.3771878 14.997027 21.483280  32.07500    1 45000
## sigma2[6897] 46.08964  9.8717144 30.675178 44.748584  69.23992    1 44560
## sigma2[6990] 29.45524  6.0321837 19.866805 28.686394  43.38169    1 45684
## sigma2[7011] 39.47482 10.5175163 23.875128 37.820573  64.68349    1 45000
## sigma2[7101] 41.85311 12.1840750 24.237028 39.730191  71.20206    1 45472
## sigma2[7172] 33.16096  7.5140550 21.569403 32.120083  50.73435    1 43921
## sigma2[7232] 46.88468  9.6019194 31.788116 45.601123  69.07229    1 44898
## sigma2[7276] 47.24516  9.6098757 32.077532 46.010865  69.51585    1 45406
## sigma2[7332] 33.46157  7.2179354 22.179040 32.487669  50.35167    1 46957
## sigma2[7341] 30.38838  6.3440396 20.468837 29.550022  45.21840    1 44577
## sigma2[7342] 42.33375  8.2046688 29.193700 41.343308  61.14980    1 44135
## sigma2[7345] 54.31211 10.7049913 37.212751 52.963423  79.02026    1 45000
## sigma2[7364] 35.57372  8.0663265 23.198048 34.423508  54.52870    1 44218
## sigma2[7635] 37.92973  7.9663986 25.450623 36.860914  56.38467    1 45202
## sigma2[7688] 21.06221  4.2749509 14.322169 20.519848  30.99096    1 45195
## sigma2[7697] 46.80730 12.6880148 28.146193 44.790827  77.15302    1 45000
## sigma2[7734] 78.51588 26.6764906 42.217812 73.596750 143.72939    1 45969
## sigma2[7890] 40.79607  8.5007120 27.355123 39.711317  60.32602    1 45000
## sigma2[7919] 48.82882 12.1798757 30.532472 46.925516  77.81224    1 44069
## sigma2[8009] 33.34019  7.1742189 22.126490 32.387149  49.96200    1 45490
## sigma2[8150] 26.70163  6.0602000 17.381521 25.843245  40.94066    1 44935
## sigma2[8165] 29.60439  6.3325262 19.665576 28.760315  44.45246    1 45925
## sigma2[8175] 39.59333 10.5725710 24.037960 37.854049  64.92911    1 45362
## sigma2[8188] 47.08454 13.2232550 27.885049 44.850434  79.22685    1 44573
## sigma2[8193] 19.68501  4.5142595 12.728704 19.026518  30.31206    1 44382
## sigma2[8202] 61.96347 15.9745772 38.233317 59.450685 100.16733    1 45825
## sigma2[8357] 48.67640 14.5584921 27.949632 46.159250  84.37082    1 43902
## sigma2[8367] 25.34600 13.0806371 10.705211 22.206306  59.02414    1 42226
## sigma2[8477] 48.94848 12.1509508 30.699126 47.132693  77.74330    1 45000
## sigma2[8531] 58.47321 13.7231198 37.521123 56.591515  90.87648    1 44572
## sigma2[8627] 43.79554  8.8865157 29.665584 42.697393  64.24942    1 45421
## sigma2[8628] 35.29145  6.6543725 24.565609 34.463447  50.50385    1 45000
## sigma2[8707] 43.13820  9.2183214 28.680666 41.952435  64.43249    1 45000
## sigma2[8775] 33.19618  7.1760534 21.911244 32.285265  49.80626    1 46404
## sigma2[8800] 44.87803 12.2748803 26.928511 42.876265  74.53022    1 46530
## sigma2[8854] 31.93908  8.8895217 18.965028 30.478818  53.38356    1 45000
## sigma2[8857] 29.73065  5.4644029 20.897571 29.105310  42.34175    1 45000
## sigma2[8874] 50.67831 12.8735408 31.455818 48.650567  81.40370    1 45122
## sigma2[8946] 43.97581  8.5301792 30.355238 42.958567  63.60680    1 45000
## sigma2[8983] 42.38942  8.8424960 28.503292 41.294313  62.75395    1 45996
## sigma2[9021] 28.43124  5.6249579 19.467669 27.698411  41.51889    1 45752
## sigma2[9104] 31.97255  6.3730835 21.822589 31.165836  46.58903    1 46148
## sigma2[9158] 41.60217  8.5227742 28.142369 40.524097  61.46057    1 45000
## sigma2[9198] 33.92014  9.5032647 20.123169 32.280653  56.71477    1 44568
## sigma2[9225] 48.51921 12.2673108 30.178061 46.697395  77.66366    1 45000
## sigma2[9292] 58.08241 21.7665828 29.372380 53.539945 112.85991    1 44575
## sigma2[9340] 51.31195 14.7098521 30.190283 48.792370  86.99259    1 43906
## sigma2[9347] 35.51519  6.9232940 24.439357 34.678935  51.41057    1 45000
## sigma2[9359] 51.88227 10.5635864 35.086644 50.541015  76.39164    1 44555
## sigma2[9397] 44.41285  9.6868769 29.379185 43.050568  66.96302    1 45000
## sigma2[9508] 44.26990 11.4188770 27.390256 42.442255  71.34067    1 45000
## sigma2[9550] 66.29829 18.9323145 38.784889 63.145955 112.07059    1 45717
## sigma2[9586] 42.60858  8.1713380 29.471864 41.628897  61.22517    1 45495
## tau00         9.08284  1.1290112  7.115918  8.999612  11.50486    1 42995

#summary(output)

การเปรียบเทียบโมเดลแบบเบส์สามารถทำได้หลายวิธีการ วิธีการหนึ่งคือการใช้ดัชนี DIC ซึ่ง JAGS คำนวณให้ DIC มีสูตรคำนวณดังนี้

กำหนดให้ deviance มีค่าเท่ากับ\(D(\theta)=-2log(p(\bf{y}|\theta))+C\) โดยที่ \(C\) คือค่าคงที่

การคำนวณค่า DIC มีหลายวิธีการ วิธีการหลักนำเสนอโดย Spiegelhalter และคณะ (2002) และ Gelman และคณะ (2004) ดังนี้

• Spiegelhalter et al. (2002) กำหนดให้ \(p_D=\overline{D(\theta})-D(\overline{\theta})\)

• Gelman และคณะ (2004) กำหนดให้ \(p_D=\frac{1}{2}\overline{Var(D(\theta))}\)

\(DIC = p_d+\overline{D(\theta)}\)

หรือ

\(DIC = D(\overline{\theta})+2p_D\)

จากสูตรในข้างต้นจะเห็นว่าโมเดลที่มีค่า DIC ต่ำกว่าจะเป็นโมเดลที่มีความเหมาะสมมากกว่า

dic.samples(fit.null, n.iter=30000, thin = 2)

## Mean deviance:  46741 
## penalty 145.5 
## Penalized deviance: 46887

dic.samples(fit.hetero, n.iter=30000, thin = 2)

## Mean deviance:  46619 
## penalty 320.9 
## Penalized deviance: 46940

package-brms

ผู้เรียนจะเห็นว่าการระบุโมเดลในภาษา JAGS ข้างต้นจำเป็นต้องอาศัยความเชี่ยวชาญพอสมควร นอกจากนี้ยังค่อนข้างใช้เวลาในการเขียนโมเดลถึงแม้จะเป็นโมเดลพื้นฐานอย่าง one-way ANOVA ในข้างต้น ปัจจุบันมี package หลายตัวที่ทำหน้าที่เป็นส่วนต่อประสานระหว่างผู้ใช้กับ engine ต่าง ๆ ซึ่งใช้ภาษาที่เรียบง่ายกว่า ทำให้การ modelling ด้วยวิธีการแบบเบส์ทำได้สะดวกและง่ายขึ้น หัวข้อนี้จะกล่าวถึง package-brms (ย่อมาจาก bayesian multilevel models using stan) ซึ่งเป็น high-level API ตัวหนึ่งของภาษา Stan รายละเอียดมีดังนี้

library(brms)
fit.brm1<-brm(formula = mathach ~ 1 + (1|schoolid), data=dat1,
              family=gaussian(),
              warmup=1000,
              iter=3000,
              chains=3)

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '0ead22ac663f58d23a13a947854bc20f' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000481 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.81 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 1: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 1: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 1: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 1: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 1: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 1: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 1: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 1: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 1: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 1: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 9.68877 seconds (Warm-up)
## Chain 1:                10.1136 seconds (Sampling)
## Chain 1:                19.8024 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '0ead22ac663f58d23a13a947854bc20f' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000322 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.22 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 2: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 2: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 2: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 2: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 2: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 2: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 2: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 2: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 2: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 2: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 9.70008 seconds (Warm-up)
## Chain 2:                10.0076 seconds (Sampling)
## Chain 2:                19.7077 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '0ead22ac663f58d23a13a947854bc20f' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000372 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.72 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 3: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 3: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 3: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 3: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 3: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 3: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 3: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 3: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 3: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 3: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 9.15719 seconds (Warm-up)
## Chain 3:                9.2879 seconds (Sampling)
## Chain 3:                18.4451 seconds (Total)
## Chain 3:

plot(fit.brm1)

summary(fit.brm1)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: mathach ~ 1 + (1 | schoolid) 
##    Data: dat1 (Number of observations: 7185) 
##   Draws: 3 chains, each with iter = 3000; warmup = 1000; thin = 1;
##          total post-warmup draws = 6000
## 
## Group-Level Effects: 
## ~schoolid (Number of levels: 160) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.95      0.19     2.60     3.33 1.00     1294     2477
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    12.64      0.25    12.17    13.15 1.00      705     1210
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     6.26      0.05     6.15     6.37 1.00    12329     4130
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

pp_check(fit.brm1, resp = "mathach")

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

bayes_R2(fit.brm1)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.1725417 0.007560215 0.1576866 0.1873903

loo(fit.brm1)

## 
## Computed from 6000 by 7185 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo -23445.7 49.4
## p_loo       146.3  2.1
## looic     46891.4 98.9
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

priors<-get_prior(mathach ~ 1 + (1|schoolid), data=dat1,
              family=gaussian())
priors

##                    prior     class      coef    group resp dpar nlpar bound
##  student_t(3, 13.1, 8.1) Intercept                                         
##     student_t(3, 0, 8.1)        sd                                         
##     student_t(3, 0, 8.1)        sd           schoolid                      
##     student_t(3, 0, 8.1)        sd Intercept schoolid                      
##     student_t(3, 0, 8.1)     sigma                                         
##        source
##       default
##       default
##  (vectorized)
##  (vectorized)
##       default

priors$prior[1]<-"normal(0,20)" #mu and sd

fit.brm2<-brm(formula = mathach ~ 1 + (1|schoolid), data=dat1,
              family=gaussian(),
              prior=priors,
              warmup=1000,
              iter=3000,
              chains=3)

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '2b3230ad39a0e3f7602f99f6d472f68b' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000464 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.64 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 1: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 1: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 1: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 1: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 1: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 1: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 1: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 1: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 1: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 1: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 8.98691 seconds (Warm-up)
## Chain 1:                9.02464 seconds (Sampling)
## Chain 1:                18.0116 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '2b3230ad39a0e3f7602f99f6d472f68b' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000285 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 2.85 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 2: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 2: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 2: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 2: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 2: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 2: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 2: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 2: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 2: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 2: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 9.12027 seconds (Warm-up)
## Chain 2:                8.97629 seconds (Sampling)
## Chain 2:                18.0966 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '2b3230ad39a0e3f7602f99f6d472f68b' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000284 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2.84 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 3000 [  0%]  (Warmup)
## Chain 3: Iteration:  300 / 3000 [ 10%]  (Warmup)
## Chain 3: Iteration:  600 / 3000 [ 20%]  (Warmup)
## Chain 3: Iteration:  900 / 3000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1001 / 3000 [ 33%]  (Sampling)
## Chain 3: Iteration: 1300 / 3000 [ 43%]  (Sampling)
## Chain 3: Iteration: 1600 / 3000 [ 53%]  (Sampling)
## Chain 3: Iteration: 1900 / 3000 [ 63%]  (Sampling)
## Chain 3: Iteration: 2200 / 3000 [ 73%]  (Sampling)
## Chain 3: Iteration: 2500 / 3000 [ 83%]  (Sampling)
## Chain 3: Iteration: 2800 / 3000 [ 93%]  (Sampling)
## Chain 3: Iteration: 3000 / 3000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 9.80854 seconds (Warm-up)
## Chain 3:                9.52466 seconds (Sampling)
## Chain 3:                19.3332 seconds (Total)
## Chain 3:

hypothesis(fit.brm1, "Intercept>0", class="sd", group="schoolid", alpha=0.05)

## Hypothesis Tests for class sd_schoolid:
##        Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (Intercept) > 0     2.95      0.19     2.65     3.27        Inf         1
##   Star
## 1    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

gamma.mcmc<-as.mcmc(as.data.frame(fit.brm1, pars="b_Intercept"))

## Warning: Argument 'pars' is deprecated. Please use 'variable' instead.

tau.mcmc<-as.mcmc(as.data.frame(fit.brm1, pars="sd")^2)

## Warning: Argument 'pars' is deprecated. Please use 'variable' instead.

sigma2.mcmc<-as.mcmc(as.data.frame(fit.brm1, variable="sigma"))^2
icc<-tau.mcmc/(tau.mcmc+sigma2.mcmc)
mean(icc)

## [1] 0.1820599

sd(icc)

## [1] 0.01892475

plotPost(icc, compVal=0.05, cenTend="median", ROPE=c(0.04,0.06), credMass=0.95)

##                  ESS      mean    median     mode hdiMass    hdiLow   hdiHigh
## Param. Val. 1323.766 0.1820599 0.1814084 0.181766    0.95 0.1445473 0.2178896
##             compVal pGtCompVal ROPElow ROPEhigh pLtROPE pInROPE pGtROPE
## Param. Val.    0.05          1    0.04     0.06       0       0       1

library(tidybayes)

## 
## Attaching package: 'tidybayes'

## The following objects are masked from 'package:brms':
## 
##     dstudent_t, pstudent_t, qstudent_t, rstudent_t

get_variables(fit.brm1)%>%head(15)

##  [1] "b_Intercept"                "sd_schoolid__Intercept"    
##  [3] "sigma"                      "r_schoolid[1224,Intercept]"
##  [5] "r_schoolid[1288,Intercept]" "r_schoolid[1296,Intercept]"
##  [7] "r_schoolid[1308,Intercept]" "r_schoolid[1317,Intercept]"
##  [9] "r_schoolid[1358,Intercept]" "r_schoolid[1374,Intercept]"
## [11] "r_schoolid[1433,Intercept]" "r_schoolid[1436,Intercept]"
## [13] "r_schoolid[1461,Intercept]" "r_schoolid[1462,Intercept]"
## [15] "r_schoolid[1477,Intercept]"

fit.brm1%>%spread_draws(r_schoolid[schoolid, Intercept]) #group-mean

## # A tibble: 960,000 × 6
## # Groups:   schoolid, Intercept [160]
##    schoolid Intercept r_schoolid .chain .iteration .draw
##       <int> <chr>          <dbl>  <int>      <int> <int>
##  1     1224 Intercept      -2.09      1          1     1
##  2     1224 Intercept      -2.74      1          2     2
##  3     1224 Intercept      -3.97      1          3     3
##  4     1224 Intercept      -1.75      1          4     4
##  5     1224 Intercept      -3.37      1          5     5
##  6     1224 Intercept      -1.55      1          6     6
##  7     1224 Intercept      -2.91      1          7     7
##  8     1224 Intercept      -3.87      1          8     8
##  9     1224 Intercept      -1.81      1          9     9
## 10     1224 Intercept      -2.17      1         10    10
## # … with 959,990 more rows

fit.brm1%>%spread_draws(b_Intercept,sd_schoolid__Intercept) #grand-mean and sqrt(tau00)

## # A tibble: 6,000 × 5
##    .chain .iteration .draw b_Intercept sd_schoolid__Intercept
##     <int>      <int> <int>       <dbl>                  <dbl>
##  1      1          1     1        12.5                   2.98
##  2      1          2     2        12.3                   3.01
##  3      1          3     3        12.5                   3.14
##  4      1          4     4        12.5                   3.00
##  5      1          5     5        12.6                   2.92
##  6      1          6     6        12.8                   2.72
##  7      1          7     7        12.5                   2.94
##  8      1          8     8        12.5                   3.15
##  9      1          9     9        12.1                   2.90
## 10      1         10    10        12.3                   2.76
## # … with 5,990 more rows

# mean and median with quantile interval
fit.brm1%>%spread_draws(b_Intercept,sd_schoolid__Intercept)%>%
              mean_qi()

## # A tibble: 1 × 9
##   b_Intercept b_Intercept.lower b_Intercept.upper sd_schoolid__Intercept
##         <dbl>             <dbl>             <dbl>                  <dbl>
## 1        12.6              12.2              13.1                   2.95
## # … with 5 more variables: sd_schoolid__Intercept.lower <dbl>,
## #   sd_schoolid__Intercept.upper <dbl>, .width <dbl>, .point <chr>,
## #   .interval <chr>

fit.brm1%>%spread_draws(b_Intercept,sd_schoolid__Intercept)%>%
              median_qi()

## # A tibble: 1 × 9
##   b_Intercept b_Intercept.lower b_Intercept.upper sd_schoolid__Intercept
##         <dbl>             <dbl>             <dbl>                  <dbl>
## 1        12.6              12.2              13.1                   2.94
## # … with 5 more variables: sd_schoolid__Intercept.lower <dbl>,
## #   sd_schoolid__Intercept.upper <dbl>, .width <dbl>, .point <chr>,
## #   .interval <chr>

random-coefficient model

priors<-get_prior(mathach~1+ses+(1+ses|schoolid), 
            data=dat, family=gaussian())
priors

##                    prior     class      coef    group resp dpar nlpar bound
##                   (flat)         b                                         
##                   (flat)         b       ses                               
##                   lkj(1)       cor                                         
##                   lkj(1)       cor           schoolid                      
##  student_t(3, 13.1, 8.1) Intercept                                         
##     student_t(3, 0, 8.1)        sd                                         
##     student_t(3, 0, 8.1)        sd           schoolid                      
##     student_t(3, 0, 8.1)        sd Intercept schoolid                      
##     student_t(3, 0, 8.1)        sd       ses schoolid                      
##     student_t(3, 0, 8.1)     sigma                                         
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default

priors$prior[1]<-"normal(0,20)" #mu and sd

fit.brm3<-brm(mathach~1+ses+(1+ses|schoolid), 
            data=dat, family="gaussian",
            prior=c(set_prior("horseshoe(1)", class="b"),
                    set_prior("student_t(1,0,5)", class="sd"),
                    set_prior("lkj(0.8)", class="cor")),
            warmup=1000, iter=6000, chains=3, thin=5,
            cores=3)

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1

## Warning: There were 323 divergent transitions after warmup. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

## Warning: The largest R-hat is 1.06, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#r-hat

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess

plot(fit.brm3)

summary(fit.brm3)

## Warning: There were 323 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See http://mc-stan.org/misc/
## warnings.html#divergent-transitions-after-warmup

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: mathach ~ 1 + ses + (1 + ses | schoolid) 
##    Data: dat (Number of observations: 7185) 
##   Draws: 3 chains, each with iter = 1200; warmup = 200; thin = 5;
##          total post-warmup draws = 600
## 
## Group-Level Effects: 
## ~schoolid (Number of levels: 160) 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)          2.21      0.15     1.94     2.53 1.00      974     2758
## sd(ses)                0.53      0.23     0.04     0.91 1.00      612     1852
## cor(Intercept,ses)    -0.16      0.29    -0.83     0.41 1.00     1929     2009
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    12.63      0.21    12.25    13.02 1.04       49       79
## ses           2.38      0.12     2.15     2.62 1.00     1471     2201
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     6.07      0.05     5.97     6.17 1.01      238      190
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

pp_check(fit.brm3, resp = "mathach")

bayes_R2(fit.brm3)

##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2206894 0.008245877 0.2042388 0.2371053

loo(fit.brm3)

## 
## Computed from 3000 by 7185 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -23238.5  51.2
## p_loo       160.3   2.5
## looic     46476.9 102.4
## ------
## Monte Carlo SE of elpd_loo is 1.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.